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APPENDIX A 
COUNTERACTING GEOMETRIC 
DISTORTIONS IN WATERMARKING 

/include <stdio.h> 
/include <stdlib.h> 
/include <math.h> 
/include <gl/gl.h> 
/include <gl/device. h> 



EXPRESS MAIL 364370875 US 



E=3 
? ? 3 



IU 

5 5= 
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FRE03 . H 



/define DC_VALUE_RATIO 0.2 

5 long number_freqs =4; 

double freqflO] = { 
1.9, 
1.654, 

10 1.476, 
1.227 

}; 

double phase [10] = { 

m 2.111, 

1.0 0.0765, 

;0 1.32, 

V\ 2.38, 
2CH 0.73, 



3.0, 
1.1 
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STAMP IT.C 



/* thia program explores the specifics of stamping an image with 'digital 
reticles' for the purpose of 
5 determining registration between a given image and the digimarc signatures, 
i.e., the registration 

requirements for finding the scale and rotation elements of the ubiquitous 
snowy image patterns */ 

10 

/include "pinecone.h" 
/include "freqs.h" 

is 

k 0 /define DISPLAY_IT 1 

| y ! /define PI 3*141592653589 

!~ /define SQRT2_2 0.707106781186 

T= /define INITIAL SCALE 10.0 

ill unsigned char *img; 
] long xdim,ydim; 

) s z Colorindex *disp; 

= u 

jfs Colorindex *scale_disp; 

?g Colorindex *pdisp; 

13 

int helplines =2; 
char *help[ J = 
{ 

30 "usage: stamp_it filename xdim ydim channels \n", 

"\n stamp_it modifies the input image to include digital reticles and 
outputs filename. stamped n 

}; 

35 

int load_scale_disp(void) { 
long i,j; 

pdisp = scale_disp+20; 
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for(j=0; j<15; j++){ 

for(i=l;i<20;i++) { 

*pdisp = (unsigned char)255 
pdisp += 20; 

} 

pdisp += 20; 

} 

pdisp » scale_disp+20+ ( 35*400 ) ; 
for(j=0; j<15; j++){ 

for(i=l;i<20;i++) { 

*pdisp = (unsigned char) 255 

pdisp += 20; 

} 

pdisp += 20; 

} 

return (0) ; 
} 

int draw_scale( double value) { 
long i; 
double dtemp; 

value *=10.0; 

memcpy (disp, scale_disp, 40000) ; 
dtemp = 20,0 * value; 
if (dtemp > 399,0)dtemp=399.0; 
pdisp = &disp[ 15*400 + (int)dtemp]; 
for(i=0;i<20;i++) { 

*pdisp ~ (unsigned char) 255; 

pdisp +=400; 

> 

rectwrite(0,0, 399, 49,disp) ; 

return (0) ; 
} 



APPENDIX A, Page 92 



WYCatfa CSTB.pa 9/2493 013) - 93 - 



double calculate_change( double distance , double scale , long which ){ 
long i; 

double value=DC_VALUE_RATIO * scale; 

for ( i=0 ; i<number_f reqs; i++ ) { 

value += scale * sin ( freq(i] * distance + phase(i) + 
(double)which * PI ); 
} 

if(value < 0.0)value = 0.0; 

return (value) ; 
} 



main( argc, argv ) 
int argc ; 

char *argv[] ; 

< 

long i, j ,go=l, button, channels, which, count; 
long temp, increment, last_middle; 

unsigned char tmp, *pimg, *pimgl , *img_out , *img_r , *img_b; 
char string[80], outfile[80]; 
FILE *inf; 

double change, current_scale = INITIAL_SCALE; 
long gid_img,gid_s tamped, gid_scale; 

if (argc I =5) { 

for( j=0; j<helplines; j++)fprintf ( stderr, n *s n , help[j]) ; 
exit( 1 ) ; 

} 



xdim « atoi(argv[2 J ) ; 

ydim = atoi ( argv( 3 J ) ; 

channels = atoi { argv [ 4 ] ) ; 

if (channels != 1 && channels 1=3) { 

fprintf( stderr, M stamp_it : channels must equal 1 for B/W 
3 for color\n" ) ; 

exit( 1 ) ; 
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} 



if (DISPLAY_IT) disp = calloc(xdira*ydim, si2eof (Color index) ) ; 

else disp = calloc(40000, sizeof (Colorindex) ) ; 
5 scale_disp = calloc(40000, sizeof (Colorindex) ) ; 

img = calloc (xdira*ydim, sizeof (unsigned char) ) ; 

img_out - calloc (xdim*ydim, sizeof (unsigned char) ) ; 

if( Idisp i! lscale_disp || I img i! ! img_out ){ 

fprintf( stderr, H stamp_it : can not allocate space\n" ) ; 
10 exit( 1 ) ; 

} 

if (channels == 3){ 

img_r - calloc ( xdim*ydim, sizeof (unsigned char) ) ; 
15*=! img_b = calloc (xdim*ydim, sizeof (unsigned char) ) ; 

if ( ! img_r j | ! img_b ) { 

fprintf( stderr, H stamp_it : can not allocate space\n p 



20= 

3 P3 



exit( 1 ) 



^ /* read in binary data into array */ 

2^ inf = fopen(argv(lJ, M r rt ); 

jj if(!inf) { 

jj, fprintf (stderr, "stamp_it: can't open %s\n" , argv{ 1 ] ) ; 

exit(l); 

} 

30 if (channels == 3){ 

f read { img_r , sizeof (unsigned char ) , xdim*ydim, inf) ; 
f read ( img, sizeof (unsigned char) ,xdim*ydim, inf) ; 
f read ( img_b, sizeof (unsigned char) , xdim*ydim, inf ) ; 

} 

35 else f read ( img, sizeof (unsigned char ), xdim*ydira, inf ) ; 

fclose(inf ) ; 

/* flip it */ 
pimg = img; 
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pimgl * &img[ xdim* (ydim-1) ] ; 
for ( i=0 ; i< ( ydim/ 2 ) ; i++ ) { 

for ( j=0; j<xdim; j++) { 

tmp = *pimg; 

*(pimg++) = *pimgl; 

*(pimgl++) = tmp; 

} 

pimgl -= (2*xdim); 

} 



if (DISPLAY_IT) { 

foreground ( ) ; 
pref size (xdim, ydim) ; 
1&, gid_img = winopen ( "Original" ) ; 

: 3 gflush(); 

!y pdisp = disp; 

sn 

; ? : pxmg = img; 

2(fi for( i=0; i< (ydim*xdim) ; i++) { 

*(pdisp++) = ( Color index) * (pimg++) ; 

j J rectwr ite (0,0, xdim-1 , ydim- 1 , d isp ) ; 

= 11 



pref size (xdim, ydim) ; 

gid_stamped = winopen( "Stamped Image. ..**); 
rectwr ite (0,0, xdim- 1 , ydim-1 , disp ) ; 



load_scale_disp( ) ; 
30 pref size(400, 50) ; 



gid_scale » winopen{ "Rough scaling..."); 



} 



35 

/* main visual feedback loop */ 
last_middle = 0; 
while(go) { 
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/* first diagonal */ 
for(i=0;i<(xdim+ydim-l) ; i++) { 
which « 0; 

/* calculate addition or subtraction to image */ 
5 change = calculate_change ( (double) i * SQRT2_2 , 

current_8cale, which ); 

if ( i < xdim ) { 

pimg = &img[ i*xdim] ; 
pimgl = &img_out [ i*xdim] ; 
10 count = i+1; 

} 

else { 

pimg = &img{xdim*ydim - ydim - xdim + i + 1]; 
pimgl = &img_out(xdim*ydim - xdim - ydim + i + 

ft 11? 

'==. count ~ xdim + ydim - i - 1; 

*5 > 

id for( j=0; j<count ; j++) { 

!fl temp = (long)*pimg + ( long) (change+0 . 5 ) ; 

20j if (temp > 255)temp = 255; 

*pimgl = (unsigned char) temp; 
pimg -= (xdim-1); 
pimgl -= (xdim-1); 

} 

} 

/* second diagonal */ 
for(i=0;i<(xdim+ydim-l) ;i++) { 
which - 1; 

/* calculate addition or subtraction to image */ 
30 change = calculate_change( (double) (xdim - 1 - i) * 

SQRT2_2 , current_scale, which ); 

if ( i < xdim ) { 

pimgl = &img_out (xdim - i - 1]; 
count = i+1; 

35 } 

else { 

pimgl = &img_out [ ( i-xdim+1) *xdim] ; 
count = xdim+ydim-i-1 ; 

} 
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for( j=0; j<count; j++) { 

temp = (long)*pimgl + ( long) (change+O. 5) 
if (temp > 2 55) temp = 255; 
*pimgl = (unsigned char) temp; 
pimgl += (xdim+1); 

} 

} 



if (DISPLAx_IT) { 

pdisp = disp; 
pimg = img_out; 
for(i=0;i<(ydim*xdim) ;i++) { 

*(pdisp++) = (Colorindex)*(pimg++) ; 

} 

winset (gid_stamped) ; 

rectwr ite ( 0, 0, xdim-1 , ydim-1 , disp) ; 

pdisp - disp; 
pimg = img; 

for(i=0;i<(ydim*xdim) ;i++) { 

* (pdisp++) = ( Color index) * (pimg++) ; 

} 

winset (gid_img) ; 

rectwr ite (0,0, xdim-1 , ydim-1 , disp) ; 



winset(gid_scale) ; 
draw_scale(current_scale) ; 



button=0; 
while( *button) { 

if( getbutton(LEFTMOUSE) ){ 

current_scale *= 1.15; 

button = 1; 

last_middle = 0; 

} 

if( getbutton(RIGHTMOUSE) ){ 

current scale *= 0.85; 
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button = 1; 
lastjniddle = 0; 

} 

if( getbutton(MIDDLEMOUSE) ){ 

while( getbutton(MIDDLEMOUSE) ); 
if( last_middle ){ 
button =1; 
go=0; 

} 

last_middle = 1; 

} 

} 

} 

/* now re-sign output image because of slight changes to master key 

/* flip output image */ 
pimg = img_out; 

pimgl - &img_out [xdim* (ydim-1) ] ; 
for ( i=0 ; i< ( ydim/2 ) ; i++ ) { 

for ( j =0 ; j <xd im ; j ++ ) { 

tmp = *pimg; 

*(pimg++) = *pimgl; 

*(pimgl++) = tmp; 

} 

pimgl -= (2*xdim); 

} 

/* write out signed image */ 
sprintf (outf ile, "%s. stamped" , argv[ 1 ] ) ; 
inf = fopen(outf ile, n w n ) ; 
if(iinf) { 

fprintf (stderr, "stamp_it: can't open %s\n" , outf ile) ; 
exit (1) ; 

} 

if (channels == 3){ 

fwrite{ img_r, sizeof (unsigned char) , xdim*ydim, inf ) ; 
fwrite( img_out , sizeof (unsigned char) , xdim*ydim, inf) ; 
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f write ( img_b, si2eof ( unsigned char ) , xdim*ydim, inf ) 

} 

else fwrite( img_out, sizeof (unsigned char ), xdim*ydim, inf ) ; 
fclose(inf ) ; 

/* free and clean up */ 
if (DISPLAY_IT) gexit(); 
f ree(img) ; 
free(img_out ) ; 
free(disp) ; 
free(scale_disp) ; 
if ( channels = 3) { 

f ree( img_r ) ; 

f reef img_b) ; 

} 
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REGISTR.C 

/* this program is a companion to stamp_it, wherein it is given a suspect 
image and it attempts to determine where 

the cross-hatch pattern resides within the suspect image, thereby determining 
the scale, roation and offset 
of the suspect image */ 

/include "pinecone.h" 
/include " f reqs . h* 

/define DISPLAY_IT 1 
/define DISP_POW 0.1 

/define MOV_AV 21 /* keep this odd please */ 

/define MAGGJTHRESHOLD 1.7 

/define ANGLE_I NCREMENT (PI/360.0) 

/define START_SCALE_ZONE 50 

/define SCALE_START 0.5 

/define SCALE_STOP 2.0 

/define SCALE_STEP 0.001 

/define MAX_BLOCK_SIZE 5 

/define PI 3.141592653589 

unsigned char *img; 

long xdira,ydim; 

Colorindex *disp; 

Colorindex *pdisp; 

long bits, f ftdim, f ft_size; 

float *ar,*ai; 

float wr [ 10000 ] ,wi[ 10000 ] ; 

int helplines =2; 
char *help{ ] = 
{ 

"usage: register filename xdim ydim channels \n", 
-\n register looks at the input image for digital reticles and 
displays registered result " 

}; 
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int shift_array( float *array,int dim) { 
int i,j; 

int dim2 = dim/2; 

int offset = dim2*dim + dim2; 

float *pl, *p2,ftmp; 

for(i=0;i<dim2;i++) { 

pi ~ &array [ i*dim) ; 

p2 = fiarray [of f set+i*dim] ; 

for( j=0; j<dim2; j++) { 

ftmp = *pl; 

*pl = *p2; 

*p2 = ftmp; 

pl++;p2++; 

} 

} 

offset = dim2*dim; 
for(i=0;i<dim2;i++) { 

pi = fiarray [dim2+i*dim] ; 
p2 = fiarray [of f set+i*dim] ; 
for ( j=0; j<dim2; j++) { 
ftmp = *pl; 
*pl = *p2; 
*p2 = ftmp; 
pl++;p2++; 

} 

} 

return(O) ; 

} 

int block_f ilter ( float *array,int dim) { 
int i,j,k,l,i2; 
float buffer [2] (4096) ; 

for(i=0;i<dim;i++)buf fer[0] [i] = array[i]; 
for(i=l;i<(dim-l) ;i++) { 
i2 = i%2; 

buf fer [ i2 ) [0] = array [ i*dim) ; 

buf fer[ i2] [dim-1] = array [ i*dim+dim-l) 
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for( j=l; j<(dim-l) ; j++){ 

buf fer[i2] ( j]=0.0; 
for(k=-l;k<2;k++) { 

for(l--l;K2;l++) { 

buffer(i2] [ j]+= array [( i+k) *dim + 

(j+D); 

} 

} 

buffer(i2] [ j]/=9.0; 

> 

i2 = (i-l)%2; 

for( j=0; j<dim; j++)array[ (i-l)*dim + j ] » buf f er [ i2 ] I j ] ; 

> 

i2 = (i-l)%2; 

for( j=0; j<dim; j++)array[ (i-l)*dim + j] = buf f er [ i2 ] { j ] ; 
return (0) ; . 

} 



/* assumes fftdim arrays ar and ai */ 
double get_mag (double x, double y){ 

long xoff,yoff; 

float *par,*pai; 

double xdist , ydist , cf , revalue, i_value , value=0 . 0 ; 

xoff = (long)x; 

yoff = (long)y; 

xdist = x - (double)xof f ; 

ydist = y - (double)yof f ; 

par » &ar[yoff * fftdim + xoff]; 
pai = &ai[yoff * fftdim + xoff]; 

cf = (1.0 - xdist) * (1.0-ydist); 
r_value = (double) * (par++) ; 
i_value = (double) * (pai++ ) ; 
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value = cf * sqrt (r_value*r_value + i_value*L_value) ; 
cf » xdist * (1.0-ydist); 
revalue » (double) *par; 
i_value 33 (double) *pai; 
5 value +» cf * sqrt (r_value*r_value + i_value*i_value) ; 

par ♦= (fftdim-1); 
pai +- (fftdim-1) ; 

cf » (1.0 - xdist) * ydist; 
10 r_value » (double)* (par++) ; 

i_value * (double) * (pai++) ; 

value += cf * sqrt (r_value*r_value + i_value*i_value) ; 
cf = xdist * ydist; 
r_value « (double) *par; 
15=| i_value = (double) *pai; 

% q value +« cf * sqrt (r_value*r_value + i_value* i_value) ; 

\Q 

return (value) ; 

a na 



= i = 

raain( argc, argv ) 

Id? 

int argc ; 

char *argv[] ; 
{ 

30 long i, j , k, go, channels, count , center; 

long dim, angle_int , total_angles, top_candidate; 

unsigned char tmp, *pimg, *pimgl, *img_out , *img_r , *img_b; 

char string[80], outfile[80); 

FILE *inf; 
35 long gid_img; 

double 

x,y, co88, sinn, angle, mag, magg( 5000 J , magg_ma[ 5000 ] , angle_vector[ 10000 ) ; 

double radius, dtmp, dscale,grey_diff, highest , f rac,highest_scale, scale; 
float *par,*pai; 
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if(argct=5) { 

for( j=0; j<helplines; j++)fprintf ( stderr, "%s", help[jj) ; 
exit( 1 ) ; 

} 

xdim = atoi(argv[2] ) ; 

ydim = atoi(argv[ 3 ] ) ; 

channels = atoi(argv[4] ) ; 

if (channels 1= 1 && channels 1=3) { 

fprintf ( stderr, "register : channels must equal 1 for B/W or 

3 for colorXn" ) ; 

exit( 1 ) ; 

} 

/* find the next power of two equal to or higher than the highest 
input dimension */ 

if (xdim > ydim) dim « xdim; 
else dim = ydim; 
fftdim = 1; go = 1; bits = 0; 
while ( go ){ 

if ( dim > fftdim ) { 
fftdim*=2; 
bits++; 

} 

else go = 0; 

} 

if (bits > 12) { 

fprintf ( stderr, "recognize : sorry, this particular program 
only accepts 4K images and lessXn" ) ; 
exit( 1 ) ; 

} 

fft_size = fftdim * fftdim; 

disp = calloc(f ft_size, sizeof (Colorindex ) ) ; 
img « calloc(xdim*ydim, sizeof (unsigned char) ) ; 
ar « calloc(f ft_size, sizeof ( float ) ) ; 
ai « calloc(f ft_size, sizeof (float ) ) ; 
if ( Idisp | j limg j j lar ! j !ai ){ 

fprintf ( stderr, "register : can not allocate space\n" ) ; 
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exit( 1 ) ; 

} 

if (channels == 3){ 

5 i- m 9_ r = calloc(xdim*ydim, sizeof (unsigned char) ) ; 

img_b = calloc (xdim*ydim, sizeof (unsigned char) ) ; 
if ( !img_r 1 | limg_b ) { 

fprintf( stderr, "register : can not allocate space\n' 

) ? 

10 exit( 1 ) ; 

} 

} 

/* read in binary data into array */ 
inf = fopen(argv[l] , n r") ; 
\Q if(iinf) { 

\Q fprintf (stderr, "register : can't open %s\n" , argv[ 1 ] ) ; 

jU exit(l); 

VA } 

2& if (channels == 3){ 

m f read ( img_r , sizeof (unsigned char) , xdim*ydim, inf) ; 

fread( img, sizeof (unsigned char ), xdim*ydim, inf ) ; 
13 f read ( img_b, sizeof (unsigned char) , xdim*ydim, inf ) ; 

!S > 

2^5 else f read ( img, sizeof (unsigned char ), xdim*ydim, inf ) ; 

■3j fclose(inf); 

/* flip it */ 
pimg = img; 
30 pimgl » &img(xdim* (ydim-1 ) ] ; 

for(i=0;i<(ydim/2) ;i++) { 

for( j=0; j<xdim; { 
tmp = *pimg; 
*(pirag++) = *pimgl; 
35 *(pimgl++) = tmp; 

} 

pimgl -= (2*xdim); 

> 
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/* copy image buffer into ar */ 
par = ar; 
pimg a img; 
for(i=*0;i<ydim;i++) { 

for ( j =0 ; j<xdim; j ++ ) { 

*(par++) = (float) *(pimg++); 

} 

par += (f ftdim-xdim) ; 

} 

/* 2d fft, in place */ 
printf ("\nforward fft... \n"); 
f ft2d(ar, ai,bits,0, wr,wi) ; 
printf ("done \n.... "); 
shift_array(ar, f ftdim) ; 
shift_array (ai, f ftdim) ; 
center = f ftdim/ 2; 

block_f ilter(ar, f ftdim) ; 
block f ilter(ai, f ftdim) ; 



if (DISPLAY_IT) { 

foreground ( ) ; 
pref size(f ftdim, f ftdim) ; 
gid_img = winopen ( "yahh" ) ; 
gf lush( ) ; 



dtmp = 

(double) ar[center*ff tdim+center+1]* (double ) ar ( center*f f tdim+center+1] 
dtmp += 

(double) ai [ center* ff tdim+center+1]* (double ) ai( center* ff tdim+center+1 J 
dscale = pow(dtmp, DISP_POW) ; 



pdisp = disp; 
par = ar; 
pai = ai; 

for ( i=0 ;i<(f ftdim* f ftdim) ;i++) { 
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dtmp = pow( (*par * *par + *pai * *pai) , 

DISP_POW ); 

dtmp *= {255.0 / dscale); 
if (dtmp>255.0)dtmp » 255.0; 
*(pdisp++) = (Colorindex) { dtmp ); 
par++;pai++; 

} 

rectwrite(0, 0, f ftdim-1, f f tdim-1 ,disp) ; 

} 



/* now search for the gross rotation axes */ 

for(angle = 0.0,angle_int = 0;angle<PI/2;angle+=ANGLE_INCREMENT, 
angie_int++) { 

coss = cos (angle); 
1§3 sinn = sin(angle); 

=.Q /* fill radial vector */ 

for(i=0;i<(f ftdim/2) ;i++) { 
radius = (double) i; 

x = (double) center - radius * coss; 
2&| y * (double) center - radius * sinn; 

!fi mag « get_mag<x, y ) ; 

5 x = (double) center + radius * sinn; 

^ y = (double) center - radius * coss; 

J J; mag += get_mag(x,y) ; 

maggfi] = mag; 

} 

/* create moving average */ 
magg_ma[MOV_AV/2] = 0.0; 
for(i=0;i<MOV_AV;i++) { 
30 magg_ma [ MOV_AV/2 J += magg [ i ] ; 

} 

magg_ma [ MOV_AV/2 ] /= ( ( double ) MOV_AV ) ; 
f or ( i= (MOV_AV/2 ) +1 ; i< ( f ftdim/2 ) - ( MOV_AV/2 ) -1 ; i++ ) { 
magg_ma [ i ] = magg_ma ( i - 1 ] ; 
35 magg_ma[i) -= ((magg[i - (MOV AV/2) - 



1 ] ) / ( double ) MOV_AV ) ; 

} 



magg_ma[i] += ( (magg [ (MOV_AV/2 ) + i ])/ (double )MOV AV) ; 
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/* within prescribed 'scale zone', calculate final number for 

this angle */ 

angle_vector[angle_int ] = 0.0; 

for { i=START_SCALE_ZONE; i< ( f f tdim/2 ) - (MOV_AV/2) -l;i++){ 
mag * magg [ i J / magg_ma [ i ] ; 
if{ mag > MAGGJTHRESHOLD ){ 

mag -= MAGG_THRESHOLD; 

angle_vector [angle_int ] (mag*mag); 

> 

} 

} 

total_angles = angle_int; 

/* sort out the TOP_CAND I DATES and find which has the best match on 
absolute scale */ 

/* chhose the highest angle_vector number for starters */ 
highest = 0.0; 

for ( angle_int=0 ; angle_int<total_angles ; angle_int++ ) { 
if ( angle_vector [ angle_int ] >highest ) { 

highest = angle_vector [ angle_int J ; 
top_candidate = angle_int; 

} 

> 

printf ( w \n\n tilt from original found = %d w , ( top_candidate/2 ) - 

45); 

coss = cos ( ANGLE_INCREMENT * (double) top_candidate) ; 
sinn = sin(ANGLE_INCREMENT * (double) top_candidate ) ; 
/* fill radial vector for this angle */ 
for(i=0;i<(f ftdim/2) ;i++) { 
radius = (double) i; 

x = (double) center - radius * coss; 
y = (double) center - radius * sinn; 
mag = get_mag(x,y); 

x = (double) center + radius * sinn; 
y = (double) center - radius * coss; 
mag += get_mag ( x , y ) ; 
magg [ i ] = mag ; 
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> 

/* create moving average */ 
raagg_ma(MOV_AV/2] = 0.0; 
for(i=0;i<MOV_AV;i++) { 

maggjna [ M0V_AV/ 2 J += magg [ i ] ; 

} 

magg_ma[MOV_AV/2] /= ( ( double )M0V_AV ); 
for ( i= (MOV_AV/2 ) +1 ; i< ( f f tdim/2 ) - (MOV_AV/2 ) -1; i++ ) { 
magg_ma [ i ] » magg_ma [ i - 1 ] ; 

magg_ma(i] -= ({magg[i - (MOV_AV/2) - 1] )/ (double )MOV_AV) ; 
magg_ma[i] += ( (magg[ (MOV_AV/2 ) + i ])/ (double )MOV_AV) ; 

} 

/* now slide the scale and find the highest point */ 
highest = 0.0; 

for(scale = SCALE_START; scale < SCALE_STOP; scale+=SCALE_STEP) { 
mag = 0.0; 

for ( j=0; j<number_f reqs; j++ ) { 

radius = scale * freq(j) * (double) fftdim / PI / 2.0 
if( (int)radius <= (l+MOV_AV/2) !! ( int ) ( radius+1) > 
((fftdim/2) - (MOV_AV/2) - 1) ); 

else { 

frac = radius - (double) ( (int) radius); 
mag += (1.0-frac) *(magg[ (int)radius ] / 

magg_ma[ (int) radius ]); 

mag += frac*(magg[ ( int ) (radius+1 ) ) / 

magg_ma[ (int) (radius+1) ]); 

} 

} 

if (mag > highest) { 

highest = mag; 
highest_scale = scale; 

} 

} 

printf ( "\n\nscale found = %lf \n n , 1 . 0/highest_scale) ; 

/* now find the exact offset and orietnation, i.e. if it is flipped, 

etc. */ 
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/* display the result at correct rotation and pixel size */ 
sleep( 1000) ; 



/* free and clean up */ 
if (DISPLAY_IT) gexit(); 
free(img) ; 
f ree(disp) ; 
10 free(ar) ; 

free(ai) ; 

if{ channels =* 3){ 

f ree ( img_r ) ; 
free ( img_b) ; 

*S } 
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